#Replication code for "Politics or prejudice: explaining individual-level hostilities in India's Hindu-Muslim conflict"
#Sebastian Schutte (sebastian@prio.org) 2019

rm(list=ls())
library(MASS)
library(stargazer)
library(coefplot)
library(xtable)
library(ggplot2)
library(flexmix)


######################
###Data preparation###
######################

#setwd('</path/to/repdata>')
data_hindu_muslim <- read.csv2("repdata_ii_final.csv",sep=";")
data_hindu_muslim_o <- data_hindu_muslim
data_hindu_muslim_o$stronger_against_outgroup <- as.ordered(data_hindu_muslim_o$stronger_against_outgroup)
data_hindu_muslim_o$violence_justified <- as.ordered(data_hindu_muslim_o$violence_justified)


data_hindu <- subset(data_hindu_muslim,religion_numeric == "hindu")
data_muslim <- subset(data_hindu_muslim,religion_numeric == "muslim")

#############################################
###Estimating models for the main analysis###
#############################################

m_prej_hindu_muslim_strato_fe <- lm(stronger_against_outgroup ~ personal_experiences_outgroup + 
                                   experienced_violence_by_outgroup + contact_with_outgroup + religion + as.factor(participant_region), data = data_hindu_muslim, weights = strato_weight)
m_secure_hindu_muslim_fe_strato <- lm(stronger_against_outgroup ~ personal_concern_outgroup + religion +  as.factor(participant_region), data = data_hindu_muslim,weights = strato_weight)

m_leaders_hindu_muslim_fe_strato <- lm(stronger_against_outgroup ~ leaders_warn_against_future_outgroup + 
                                  leaders_warn_against_past_outgroup + religion + as.factor(participant_region), data = data_hindu_muslim, weights = strato_weight)

m_full_hindu_muslim <- lm(stronger_against_outgroup ~ age + gender + years_in_school + 
                            family_relocate + religious_identification + 
                            #Political controls
                            people_wrong_should_stay_out +  democracy_important + ingroup_shares_view + outgroup_creates_problems +
                            #Security
                            personal_concern_outgroup + 
                            #Elite
                            leaders_warn_against_future_outgroup + leaders_warn_against_past_outgroup + 
                            #Prejudice
                            personal_experiences_outgroup + experienced_violence_by_outgroup + 
                            contact_with_outgroup + city_riots_pre + religion #+ ab_filter, data = data_hindu_muslim)
                          , data = data_hindu_muslim)

m_full_hindu_muslim_fe <- lm(stronger_against_outgroup ~ age + gender + years_in_school + 
                            family_relocate + religious_identification + 
                            #Political controls
                            people_wrong_should_stay_out +  democracy_important + ingroup_shares_view + outgroup_creates_problems +
                            #Security
                            personal_concern_outgroup + 
                            #Elite
                            leaders_warn_against_future_outgroup + leaders_warn_against_past_outgroup + 
                            #Prejudice
                            personal_experiences_outgroup + experienced_violence_by_outgroup + 
                            contact_with_outgroup + city_riots_pre + religion  +  as.factor(participant_region), data = data_hindu_muslim)

m_full_hindu_muslim_fe_strato <- lm(stronger_against_outgroup ~ age + gender + years_in_school + 
                            family_relocate + religious_identification + 
                            #Political controls
                            people_wrong_should_stay_out +  democracy_important + ingroup_shares_view + outgroup_creates_problems +
                            #Security
                            personal_concern_outgroup + 
                            #Elite
                            leaders_warn_against_future_outgroup + leaders_warn_against_past_outgroup + 
                            #Prejudice
                            personal_experiences_outgroup + experienced_violence_by_outgroup + 
                            contact_with_outgroup + city_riots_pre + religion  +  as.factor(participant_region), data = data_hindu_muslim, weights = strato_weight)


m_vio_hindu_muslim_full <- lm(violence_justified ~ stronger_against_outgroup + age + gender + years_in_school + 
                              family_relocate + religious_identification + outgroup_creates_problems + leaders_warn_against_future_outgroup +
                              leaders_warn_against_past_outgroup + personal_concern_outgroup + personal_experiences_outgroup + 
                              experienced_violence_by_outgroup + contact_with_outgroup + people_wrong_should_stay_out + 
                              democracy_important + ingroup_shares_view + city_riots_pre + religion, 
                              data = data_hindu_muslim)

#Main regression table
stargazer(m_full_hindu_muslim, m_full_hindu_muslim_fe, m_full_hindu_muslim_fe_strato, m_secure_hindu_muslim_fe_strato, 
                m_leaders_hindu_muslim_fe_strato, m_prej_hindu_muslim_strato_fe, m_vio_hindu_muslim_full,omit = c("as.factor"),type = "text")


#########################################################
###Coefficient plot for the main explanatory variables###
#########################################################

# Put model estimates into temporary data.frames:
model1Frame <- data.frame(Variable = rownames(summary(m_full_hindu_muslim)$coef),
                          Coefficient = summary(m_full_hindu_muslim)$coef[, 1],
                          SE = summary(m_full_hindu_muslim)$coef[, 2],
                          Model = "Basic model")
#Get the number of predictors:
k_vars <- length(m_full_hindu_muslim$coefficients)
model2Frame <- data.frame(Variable = rownames(summary(m_full_hindu_muslim_fe)$coef)[1:k_vars],
                          Coefficient = summary(m_full_hindu_muslim_fe)$coef[1:k_vars, 1],
                          SE = summary(m_full_hindu_muslim_fe)$coef[1:k_vars, 2],
                          Model = "With Regional FE")
model3Frame <- data.frame(Variable = rownames(summary(m_full_hindu_muslim_fe_strato)$coef)[1:k_vars],
                          Coefficient = summary(m_full_hindu_muslim_fe_strato)$coef[1:k_vars, 1],
                          SE = summary(m_full_hindu_muslim_fe_strato)$coef[1:k_vars, 2],
                          Model = "With Regional FE, Post-stratified")

# Combine these data.frames
allModelFrame <- data.frame(rbind(model1Frame, model2Frame, model3Frame))

allModelFrame[,1] <- rownames(summary(m_full_hindu_muslim)$coef)

#Only retain variables of interest for the plot
allModelFrame <- allModelFrame[which(grepl("concern",allModelFrame[,1]) |
                                       grepl("future_",allModelFrame[,1]) |
                                       grepl("person",allModelFrame[,1]) | 
                                       grepl("creat",allModelFrame[,1]) |
                                       grepl("share",allModelFrame[,1]) | 
                                       grepl("reloc",allModelFrame[,1])),]

allModelFrame$Variable[grepl("family_relocate",allModelFrame$Variable)] <- "C: Family relocation"
allModelFrame$Variable[grepl("concern",allModelFrame$Variable)] <- "H1: Personal concerns"
allModelFrame$Variable[grepl("create",allModelFrame$Variable)] <- "C: Out-group \"creates problems\""
allModelFrame$Variable[grepl("future_",allModelFrame$Variable)] <- "H2: Leaders point to future problems"
allModelFrame$Variable[grepl("shares_view",allModelFrame$Variable)] <- "C: In-group shares view"
allModelFrame$Variable[grepl("person",allModelFrame$Variable)] <- "H3: Personal experiences with out-group"
rownames(allModelFrame) <- c(1:nrow(allModelFrame))

#Center plot title
theme_update(plot.title = element_text(hjust = 0.5))

# Specify the width of your confidence intervals
interval1 <- -qnorm((1-0.9)/2)  # 90% multiplier
interval2 <- -qnorm((1-0.95)/2)  # 95% multiplier

# Plot
zp1 <- ggplot(allModelFrame, aes(shape=Model))
zp1 <- zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
zp1 <- zp1 + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
                                ymax = Coefficient + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                 ymax = Coefficient + SE*interval2),
                             lwd = 1/2, position = position_dodge(width = 1/2),
                             fill = "WHITE")
zp1 <- zp1 + coord_flip() + theme_bw()
zp1 <- zp1
print(zp1)


####################################
###Estimate finite mixture models###
####################################

#religion dummy
data_fmm <- as.data.frame(c(rep(1,nrow(data_hindu)),rep(0,nrow(data_muslim))))
#dp
data_fmm <- cbind(data_fmm,c(data_hindu$hindus_stronger_against_muslims,data_muslim$muslims_stronger_against_hindus))
#concerns
data_fmm <- cbind(data_fmm,c(data_hindu$personal_concern_muslims,data_muslim$personal_concern_hindus))
#leaders_future
data_fmm <- cbind(data_fmm,c(data_hindu$leaders_warn_against_future_muslims,data_muslim$leaders_warn_against_future_hindus))
#leaders_past
data_fmm <- cbind(data_fmm,c(data_hindu$leaders_warn_against_past_muslims,data_muslim$leaders_warn_against_past_hindus))
#experiences
data_fmm <- cbind(data_fmm,c(data_hindu$personal_experiences_muslims,data_muslim$personal_experiences_hindus))
#violence
data_fmm <- cbind(data_fmm,c(data_hindu$experienced_violence_by_muslims,data_muslim$experienced_violence_by_hindus))
#contact
data_fmm<- cbind(data_fmm,c(data_hindu$contact_with_muslims,data_muslim$contact_with_hindus))
#others_share_view
data_fmm<- cbind(data_fmm,c(data_hindu$hindus_share_view,data_muslim$muslims_share_view))
#outgroup_problem
data_fmm<- cbind(data_fmm,c(data_hindu$muslims_create_problems,data_muslim$hindus_create_problems))
#location information
data_fmm<- cbind(data_fmm,c(data_hindu$participant_latitude,data_muslim$participant_latitude))
data_fmm<- cbind(data_fmm,c(data_hindu$participant_longitude,data_muslim$participant_longitude))

names(data_fmm) <- c("hindu","stronger_stance","concerns","leaders_future","leaders_past","experiences",
                     "exp_violence","contact","others_share_view","outgroup_causes_probs","latitude","longitude")

#Listwise deletion of NAs 
data_fmm <- na.omit(data_fmm)


m_sec <- flexmix(stronger_stance ~  hindu + concerns, 
                  data = data_fmm, model = FLXMRglm(family = "gaussian"), k = 3)
m_lead <- flexmix(stronger_stance ~  hindu + leaders_future + leaders_past, 
                 data = data_fmm, model = FLXMRglm(family = "gaussian"), k = 3)
m_prej <- flexmix(stronger_stance ~  hindu +  experiences + exp_violence + contact,
                 data = data_fmm, model = FLXMRglm(family = "gaussian"), k = 3)

#Setting random seed, increasing max iterations
set.seed(123)
mycont <- list(iter.max = 1000, tol = 0.001, class = "r")

model_fmm <- FLXMRglmfix(nested=list(k=c(1,1,1), formula=c(~ hindu + concerns, 
                                                        ~ hindu + leaders_future + leaders_past, 
                                                        ~ hindu +  experiences + exp_violence + contact)))

result_fmm <- stepFlexmix(stronger_stance ~ 1, k=3, model=model_fmm, data=data_fmm, nrep=500, control = as(mycont, "FLXcontrol"))

sec_obs  <- data_fmm[which(result_fmm@cluster == 1),]
lead_obs <- data_fmm[which(result_fmm@cluster == 2),]
prej_obs <- data_fmm[which(result_fmm@cluster == 3),]
cat(nrow(sec_obs),nrow(lead_obs),nrow(prej_obs),"\n")

par(mfrow=c(3,1))
plot(hist(sec_obs$stronger_stance, breaks =  max(sec_obs$stronger_stance) - min(sec_obs$stronger_stance), plot = F), 
     col=rgb(0.8,0.8,0.8,0.5), xlim=c(1,5), main = "Support for stronger stance for \n observations from the \"security dillemma\" model", xlab="")  # sec
plot(hist(lead_obs$stronger_stance, breaks =  max(lead_obs$stronger_stance) - min(lead_obs$stronger_stance), plot = F), 
     col=rgb(0.4,0.4,0.4,0.5), xlim=c(1,5), main = "Support for stronger stance for \n observations from the \"elite manipulation\" model", xlab="")  # lead  
plot(hist(prej_obs$stronger_stance, breaks =  max(prej_obs$stronger_stance) - min(prej_obs$stronger_stance), plot = F), 
     col=rgb(0.1,0.1,0.1,0.5), xlim=c(1,5), main = "Support for stronger stance for \n observations from the \"prejudice\" model", xlab="")#"Agreement to a stronger stance")  # prej

############################
###Supplementary material###
############################


##############
#Descriptives#
##############
matrix_vars_hindus <- c("age","gender","years_in_school","family_relocate","religious_identification",
                        "people_wrong_should_stay_out_h","democracy_important_h","hindus_share_view","muslims_create_problems",
                        "personal_concern_muslims","leaders_warn_against_future_muslims","leaders_warn_against_past_muslims",
                        "personal_experiences_muslims","experienced_violence_by_muslims","contact_with_muslims")

matrix_vars <- matrix_vars_hindus
corr_mat <- matrix(nrow=length(matrix_vars),ncol=length(matrix_vars))
colnames(corr_mat) <- matrix_vars
rownames(corr_mat) <- matrix_vars
for (i in 1:length(matrix_vars)){
  cvar <-  colnames(corr_mat)[i]
  cvals <- data_hindu[,names(data_hindu) == cvar]
  for(j in 1:i){
    rvar <-  rownames(corr_mat)[j]
    rvals <- data_hindu[,names(data_hindu) == rvar] 
    corr_mat[i,j] <- cor(cvals,rvals,use = "complete.obs")
    
  }
}
xtable(round(corr_mat,3))

matrix_vars_muslims <- c("age","gender","years_in_school","family_relocate","religious_identification",
                         "people_wrong_should_stay_out_m","democracy_important_m","muslims_share_view","hindus_create_problems",
                         "personal_concern_hindus","leaders_warn_against_future_hindus","leaders_warn_against_past_hindus",
                         "personal_experiences_hindus","experienced_violence_by_hindus","contact_with_hindus")

matrix_vars <- matrix_vars_muslims
corr_mat <- matrix(nrow=length(matrix_vars),ncol=length(matrix_vars))
colnames(corr_mat) <- matrix_vars
rownames(corr_mat) <- matrix_vars
for (i in 1:length(matrix_vars)){
  cvar <-  colnames(corr_mat)[i]
  cvals <- data_muslim[,names(data_muslim) == cvar]
  for(j in 1:i){
    rvar <-  rownames(corr_mat)[j]
    rvals <- data_muslim[,names(data_muslim) == rvar] 
    corr_mat[i,j] <- cor(cvals,rvals,use = "complete.obs")
    
  }
}
xtable(round(corr_mat,3))

dist_vars <- c(matrix_vars_hindus,"city_riots_pre")
col_names <- c("Variable","Mean","Standard deviation","Minimum","Maximum","Observations","Definition")
dist_mat <- matrix(nrow=length(dist_vars),ncol=length(col_names))
rownames(dist_mat) <- dist_vars
colnames(dist_mat) <- col_names
for(i in 1:length(dist_vars)){
  dist_mat[i,2] <- mean(data_hindu[,names(data_hindu) == dist_vars[i]],na.rm=T)	
  dist_mat[i,3] <- sd(data_hindu[,names(data_hindu) == dist_vars[i]],na.rm=T)
  dist_mat[i,4] <- min(data_hindu[,names(data_hindu) == dist_vars[i]],na.rm=T)
  dist_mat[i,5] <- max(data_hindu[,names(data_hindu) == dist_vars[i]],na.rm=T)
}
xtable(round(dist_mat,3))


dist_vars <- c(matrix_vars_muslims,"city_riots_pre")
col_names <- c("Variable","Mean","Standard deviation","Minimum","Maximum","Observations","Definition")
dist_mat <- matrix(nrow=length(dist_vars),ncol=length(col_names))
rownames(dist_mat) <- dist_vars
colnames(dist_mat) <- col_names
for(i in 1:length(dist_vars)){
  dist_mat[i,2] <- mean(data_muslim[,names(data_muslim) == dist_vars[i]],na.rm=T)	
  dist_mat[i,3] <- sd(data_muslim[,names(data_muslim) == dist_vars[i]],na.rm=T)
  dist_mat[i,4] <- min(data_muslim[,names(data_muslim) == dist_vars[i]],na.rm=T)
  dist_mat[i,5] <- max(data_muslim[,names(data_muslim) == dist_vars[i]],na.rm=T)
}
xtable(round(dist_mat,3))

####################
#Interview duration#
####################
par(mfrow=c(1,1))
plot(density(subset(data_hindu_muslim,religion < 3)$duration/60,na.rm=T),
     main="Survey completion in minutes\n for Hindus and Muslims")

#######################################
#Robustness check excluding Tamil Nadu#
#######################################

m_full_hindu_muslim_tn <- lm(stronger_against_outgroup ~ age + gender + years_in_school + 
                               family_relocate + religious_identification + 
                               #Political controls
                               people_wrong_should_stay_out +  democracy_important + ingroup_shares_view + outgroup_creates_problems +
                               #Security
                               personal_concern_outgroup + 
                               #Elite
                               leaders_warn_against_future_outgroup + leaders_warn_against_past_outgroup + 
                               #Prejudice
                               personal_experiences_outgroup + experienced_violence_by_outgroup + 
                               contact_with_outgroup + city_riots_pre + religion, data = subset(data_hindu_muslim,participant_region != 'Tamil Nadu'))

stargazer(m_full_hindu_muslim_tn,type = "text")


###############################################
#Ordered logisitc regression for main analysis#
###############################################
mo_prej_hindu_muslim <- polr(stronger_against_outgroup ~ personal_experiences_outgroup + 
                                      experienced_violence_by_outgroup + contact_with_outgroup + religion, data = data_hindu_muslim_o)
mo_secure_hindu_muslim <- polr(stronger_against_outgroup ~ personal_concern_outgroup + religion, data = data_hindu_muslim_o)
mo_leaders_hindu_muslim <- polr(stronger_against_outgroup ~ leaders_warn_against_future_outgroup + 
                                         leaders_warn_against_past_outgroup + religion, data = data_hindu_muslim_o)

mo_full_hindu_muslim <- polr(stronger_against_outgroup ~ age + gender + years_in_school + 
                            family_relocate + religious_identification + 
                            #Political controls
                            people_wrong_should_stay_out +  democracy_important + ingroup_shares_view + outgroup_creates_problems +
                            #Security
                            personal_concern_outgroup + 
                            #Elite
                            leaders_warn_against_future_outgroup + leaders_warn_against_past_outgroup + 
                            #Prejudice
                            personal_experiences_outgroup + experienced_violence_by_outgroup + 
                            contact_with_outgroup + city_riots_pre + religion #+ ab_filter, data = data_hindu_muslim)
                          , data = data_hindu_muslim_o)


mo_vio_hindu_muslim_full <- polr(violence_justified ~ stronger_against_outgroup + age + gender + years_in_school + 
                                family_relocate + religious_identification + outgroup_creates_problems + leaders_warn_against_future_outgroup +
                                leaders_warn_against_past_outgroup + personal_concern_outgroup + personal_experiences_outgroup + 
                                experienced_violence_by_outgroup + contact_with_outgroup + people_wrong_should_stay_out + 
                                democracy_important + ingroup_shares_view + city_riots_pre + religion, 
                              data = data_hindu_muslim_o)

#Main regression replicated with ordered logistic regression models
stargazer(mo_full_hindu_muslim, mo_secure_hindu_muslim, mo_leaders_hindu_muslim, mo_prej_hindu_muslim, mo_vio_hindu_muslim_full,omit = c("as.factor"),type = "text")

